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ABSTRACT 

We study three processes that shape the nuclear composition of GRB fireballs: (1) 
neutronization in the central engine, (2) nucleosynthesis in the fireball as it expands 
and cools, and (3) spallation of nuclei in subsequent internal shocks. The fireballs are 
found to have a neutron excess and a marginally successful nucleosynthesis. They are 
composed of free nucleons, a-particles, and deuterium. A robust result is the survival of 
a significant neutron component, which has important implications. First, as shown in 
previous works, neutrons can lead to observable multi-GeV neutrino emission. Second, 
as we show in an accompanying paper, neutrons impact the explosion dynamics at radii 
up to 10 17 cm and change the mechanism of the GRB afterglow emission. 

Subject headings: accretion, accretion disks — cosmology: miscellaneous — dense 
matter — gamma rays: bursts — nuclear reactions, nucleosynthesis 

1. Introduction 

Cosmological gamma-ray bursts (GRBs) are powerful explosions in distant galaxies. A physical 
picture of these phenomena has emerged in the past decade (see Meszaros 2002 for a review) : the 
GRBs are generated by compact, dense, and energetic engines, and they are likely related to black 
hole formations. The typical mass of the central engine is believed to be a few M@, its size is 
10 6 — 10 7 cm, and its temperature is kT =1 — 10 MeV. The engine ejects a hot outflow ("fireball") 
made of free nucleons, e ± pairs, trapped blackbody radiation, and magnetic fields. 

The initial nuclear composition of the fireball and its evolution during expansion turn out 
crucial for the explosion physics. In particular, neutrons were shown to be an important component 
in the explosion (Derishev, Kocharovsky, & Kocharovsky 1999) and an extreme case of fireballs 
with a large neutron excess was studied (Fuller, Pruet, & Abazajian 2000). Particular attention 
was given to the relative motion of the neutron and ion components of the fireball, which can lead 
to observable multi-GeV neutrino emission (Bahcall & Meszaros 2000; Meszaros &: Rees 2000). 
The two components can be ejected with substantially different velocities, and this can affect the 
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observed afterglow emission of GRBs (Pruet & Dalai 2002; Bulik, Sikora, k, Moderski 2002). In an 
accompanying paper (Beloborodov 2002, hereafter Paper 2) we show that the presence of neutrons 
crucially changes the fireball interaction with an external medium and implies a new mechanism of 
the afterglow production. 

The present paper focuses on processes that shape the nuclear composition of GRB fireballs. 
In § 2 we study the production of neutrons by the central engine and the resulting neutron-to- 
proton ratio in the ejected material (scale R ~ tq 55 10 7 cm). In § 3, we calculate nuclear reactions 
in the expanding fireball and find abundances of survived free nucleons and synthesized helium 
(R = 10 7 — 10 9 cm). Similar nucleosynthesis calculations, with different codes, have been done 
recently by Lemoine (2002) and Pruet, Guiles, & Fuller (2002). In § 4, we study spallation of helium 
at later stages when internal motions develop in the fireball (i? = 10 9 — 10 12 cm). The subsequent 
dynamics of GRB blast waves with the survived neutron component (R = 10 15 — 10 17 cm) is 
investigated in Paper 2. 



There are various models for the central engines of GRBs. It may be a young neutron star 
whose rotational energy is emitted in a magnetized wind, a neutron star merger, or a massive star 
collapse. The latter two scenarios proceed via the formation of a black hole of mass M ~ M & 
and subsequent disk-like accretion of a comparable mass. The baryonic component of the ejected 
fireball is then picked up from the accretion disk. 

The central engines are sufficiently dense for the electron capture reaction. We will calculate 
at what densities and temperatures this process creates a neutron excess (neutron-to-proton ratio 
above unity), and then show that GRB engines are likely to satisfy these conditions. We will 
illustrate with the accretion-type models of GRBs, where the matter density can be relatively low 
and neutronization is most questionable. 



Consider a dense, p > 10 7 g cm~ 3 , and hot, kT > m e c 2 , matter. The rates of photon emission 
and absorption are huge and the matter is filled with Planckian radiation. Also, the rates of pair 
creation and annihilation (7 + 7 <-> e~ + e + ) are huge, and the pairs are in perfect thermodynamic 
equilibrium with the baryonic matter and radiation. The number densities are 



Here p is particle momentum in units of m e c and f± is the Fermi-Dirac occupation function, 
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Neutronization 



2.1. The equilibrium Y e 
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where 9 = kT/m e (? and fi = /x_ = —p + is the electron chemical potential in units of m e c 2 
(thermodynamic equilibrium of with radiation implies p++fi- = 0). In addition to equation (1), 
we have the charge-neutrality condition, 

n_-n + = Y e — , (3) 
m p 

where Y e is the proton-to-nucleon ratio, which would equal the electron-to-nucleon ratio in the 
absence of e^ pairs. Equations (1) and (3) determine p and n± for given T, p, and Y e . The 
electrons become degenerate if p exceeds 9, which happens below the characteristic degeneracy 
temperature, 

»- / \ 1/3 / \ 1/3 

** = ^(£) • kT ^ = 7J {w^) MeV - (4) 

Degeneracy exponentially suppresses the positron density, n + /n_ ~ exp(—p/9), because then 
are created only at energies E/m e c 2 ~ p > 9, in the exponential tail of the thermal distribution. 

At temperatures and densities under consideration, the baryonic matter is in nuclear statistical 
equilibrium, and it is dominated by free nucleons in the unshadowed region of Figure 1. The 
boundary of this region has been calculated with the Lattimer-Swesty Equation of State code 
(Lattimer & Swesty 1991). The free protons and neutrons can capture e~ and e + via charged 
current reactions, 

e^+p^n + v, e + +n^p J r v. (5) 

These reactions can rapidly convert protons into neutrons and neutrons back into protons, and 
establish an equilibrium Y e = n p /(n n + n p ), where n p and n n are number densities of protons and 
neutrons, respectively. 4 We now calculate the equilibrium Y e (T,p), and in § 2.2 we will show how 
it applies to GRB central engines. 

The exact equilibrium Y e depends on whether the opposite reactions — reabsorption of the 
emitted v and v — are also significant. We first consider the ^-transparent case, where reabsorption 
can be neglected, and then address the i/-opaque case. 



2.1.1. Neutrino-transparent matter 

The rates of e~ and e + capture can be derived from the standard electro- weak theory (e.g., 
Shapiro & Teukolsky 1983, Bruenn 1985). We assume the nucleons to be non-degenerate, and then 
the rates are 

/>oo 

n e - p = Kn p / /_ (oj + Q) (w + Q) 2 
Jo 



4 The neutron decay n^p + e +i/isa slow process on GRB timescales and it is neglected here. 
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n = Kn n f+(u - Q)(u - Qf 
JQ+l 
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(o;-Q) 2 . 



(7) 



where to is neutrino energy in units of m e c 2 , Q = (m n — m p )/m e = 2.531, and K ps 6.5 x 10 4 s 1 . 
The constant K can be expressed in terms of the mean lifetime of neutrons with respect to /3-decay, 
Tp srf 900 s, as K « (l^)- 1 (Shapiro & Teukolsky 1983, pp. 316, 524). 

An equilibrium Y" e is established when the rates of e~ and e + captures are equal, 5 



n 



e~p- 



(8) 



Equations (1), (3), and (8) determine Y e for given p and T. Contours of the function Y e (T,p) on 
the T — p plane are shown in Figure 1. 

The transition from a proton excess (Y e > 0.5) to a neutron excess (Y e < 0.5) takes place in 
the region of mild degeneracy where p < 6. We now focus on this region and derive the equilibrium 
Y e analytically. At p < 6 and 6 > Q + 1 equations (6) and (7) simplify, 



n, 



e~p 
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(10) 



where ("(5) = 1.037 is Riemann ^-function. We neglected here next-order terms 0(Q 2 /8 2 ), 0(p 2 /6 2 ), 
and 0[(Q + l) 5 /6> 5 ], and used the formula f™(e x + l)~ 1 x n dx = (1 - 2~ n )r(n + l)C(n + 1) with 
T(n + 1) = n\ for integer n. Equating the two rates, we get the equilibrium Y e = n p /(n n + n p ), 



Y P 



1 



+ 



7vr 4 (Q/2 - p) 
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1 + 0.487 ig^zjO. 

2 p 



(11) 



In the non-degenerate limit, p/0 ^ 0, this gives Y" e = 0.5 + 0.616/0 > 0.5 and implies a proton 
excess, which is due to the positive difference Q between the neutron and proton masses. A very 
mild degeneracy p = Q/2 < 9 is sufficient to drive Y e below 0.5. This happens because the e + 
density is reduced by the degeneracy effects and the e~ capture becomes preferential. 



It is instructive to write the e and e + densities using the linear expansion of equation (1) in 



p/9, 



n± 



7T 2 A 3 



:C(3)0 3 



^ 2 



p < 



(12) 



5 In the ^-transparent regime, neutrinos arc sometimes prescribed a zero chemical potential, and the balance 
fi + Hp = fin + [i v is used with fi v = to determine the equilibrium Y e . In fact, the balance of chemical potentials 
does not hold because the neutrinos are out of thermodynamic equilibrium. This balance would be valid only in the 
cold limit T -> (Landau & Lifshitz 1980). 
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kT, MeV 



Fig. 1. — Contours of the equilibrium Y e (T, p) on the T — p plane for ^-transparent matter. Dashed 
line shows the degeneracy temperature (eq. 4). In the shadowed region, the baryonic matter is 
dominated by composite nuclei, and the calculations based on rates (6) and (7) are not valid. 
The n/p-ratio equals (1 — Y e )/Y e , and the neutrons dominate over protons where Y e < 0.5. The 
analytically calculated boundary Y e = 0.5 (eq. 14) is also plotted here by a dotted line, which 
perfectly coincides with the numerically found contour Y e = 0.5. 
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where A = h/m e c = 3.862 x 10~ n cm and C(3) = 1.202. Then equation (3) yields 



" = 3F ^- < 13 > 



The condition p > Q/2 that defines the neutron-excess region on the T — p plane (Y e < 0.5) can 
now be written as 9 < 9 n (p), where 

1 /2 

9 n = , kT n = 33p}{ 2 MeV. (14) 

\ yip / 

This simple formula perfectly coincides with the numerical results (Fig. 1). That a very mild 
electron degeneracy is sufficient to create a neutron excess can also be seen by comparing T n with 

4.3^°. (15) 



T n _ , .,„J (i 

Tdeg 



A useful explicit formula for the equilibrium Y e (T,p) in ^-transparent matter with mild de- 
generacy is derived from equations (11) and (13), 

It agrees with the numerical calculations shown in Figure 1 with a high accuracy, 5Y e /Y e < 1% at 
Y e > 0.35. In a more degenerate region, where the equilibrium Y e < 0.35, the formula can still be 
used, though its error increases to 10% at Y e = 0.2 and 30% at Y e = 0.1. 



2.1.2. Neutrino- opaque matter 

If the matter is opaque to the emitted neutrinos, a complete thermodynamic equilibrium is 
established. A detailed balance now holds, e~ + p <-> n + u and e + + n <-> p + V, and the equilibrium 
Y e is determined by the condition 

p + p p = p n + p u , (17) 

where p, p p , p n , and p u are chemical potentials (in units of m e c 2 ) of the electrons, protons, neutrons, 
and neutrinos, respectively, all including the particle rest-mass energy. The antineutrinos have 
chemical potential p v = —p u , so that p+ + p n = p p + p v is also satisfied. The neutrons and protons 
have Maxwellian distributions, which gives n n n~ l = ex.p[(p n m e c 2 — m n c 2 )/kT] exp[(— p p m e c 2 + 
m p c 2 )/kT] and 

-Mp = 9\n{n n /n p ) + Q. (18) 

The thermalized v and v obey Fermi-Dirac statistics, and they are described in the same way 
as e ± (see eqs. [1] and [2]); the only difference is that the statistical weight of energy states is one 
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for neutrinos and two for electrons. The neutrino chemical potential vanishes if v and v have equal 
densities, n v = rip. If, however, the matter emits non-equal numbers of v and v (its Y e is changing) 
then n u / rip and p u / 0. 

Suppose N v neutrinos and Np anti-neutrinos have been emitted per nucleon. This causes a 
change of Y e , 

Y e - Y e ° = Np- N v , (19) 

where Y" e ° is an initial value that the matter had before the neutrino emission. If all the emitted 
neutrinos remain trapped in the matter then (n u — rip) /rib = Y e — Y e °, where rib = n n + n p is 
the total nucleon density. If a fraction x of the emitted neutrino diffused out of the matter then 
(n u - rip)/n b = (1 - x)(Y e - Y" e °) and 

Y° 



Ti- — n + v J \ Y e 



= (!-*) 1-^- ■ (20) 



Let us first consider the case x — > 1 (efficient neutrino cooling). Then — rip| <C n_ — n + and the 
neutrino chemical potential can be neglected compared to that of the electrons. Thus, the chemical 
equilibrium in ^-cooled, ^-opaque matter reads 

H = e\og(n n /n p ) + Q. (21) 

Equation (21) combined with equations (1) and (3) determines an equilibrium Y e as a function of 
T and p. Note however that equation (21) assumes free nucleons, which is invalid in the shadowed 
region of the T — p plane shown in Figure 1 (the boundary of this region is also plotted in Fig. 2). 
The chemical balance can be extended to this region if p n and p p are corrected for the heavy 
nuclei formation, which we do using the Lattimer-Swesty Equation of State code. Then we find the 
equilibrium Y e that satisfies p = p n — p p . The results are shown in Figure 2. 

The neutron excess boundary, Y e = 0.5, lies in the region where all nucleons are free and the 
electrons are mildly degenerate, p < 9. Equation (21) applies here and shows that this boundary is 
defined by p = Q. Using the linear expansion in p,/6 (eq. 13), we get a simple equation for Y e (9, p), 



■ e 



In particular, one sees that Y e = 0.5 corresponds to 



o \3 \ 1/2 

P ' , kT n = 23.1p\{ 2 MeV. (23) 



2Qm p 



It coincides exactly with the contour Y e = 0.5 calculated with the Lattimer-Swesty code (Fig. 2). 
Note also that T n = 3T dcg p\{ 6 . 

Finally, we address the regime where the neutrinos are not only thermalized but also trapped 
in the matter (x « in eq. [20]). It can happen in GRB accretion flows with high accretion rates, 
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kT, MeV 



Fig. 2. — Contours of the equilibrium Y e (T, p) on the T — p plane for i/-opaque matter with p u = 0. 
Dashed line shows the degeneracy temperature (eq. 4). Thick solid curve is the boundary of the 
free-nucleon region (same as in Fig. 1). Above this curve, the composite nuclei dominate, which 
prefer equal numbers of neutrons and protons, and therefore the contours Y e = const bend upward. 
The analytically calculated boundary Y e = 0.5 (eq. 23) is shown by dotted line; it coincides with 
the numerically found contour Y e = 0.5. 
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M > 1Mq/s (e.g. Di Matteo, Perna, & Narayan 2002). Then //„ should not be neglected in the 
chemical balance. Like equation (13) we derive for neutrinos 

6A 3 (nj, - z/ p ) 



Mi/ 



e 2 



(24) 



and substitute n v — v v = (Y e — Y®)(p/m p ). The chemical balance (i — fi v = fi n — fi p now reads 



3(2y e ° - Y e ) pA 3 



2 



6 log 



l-Y e 



Y P 



+ Q, 



and Y e = 0.5 corresponds to 



3(2y e ° - 0.5)A 3 /9 



1/2 



(25) 



(26) 



At Y e ° = 0.5 it coincides with equation (23) as it should: Y e = 0.5 = Y® requires that v and v are 
emitted in equal numbers and \i v = 0. 



2.1.3. Effects of magnetic fields 

Magnetic fields have been neglected in the above calculations, which is a valid approximation as 
long as the field does not affect the particle distribution functions. Magnetic fields can be generated 
in the GRB central engines by dynamo, and their energy is likely a fraction e# < 1 of the total 
energy density w, so that 

B = ^8tt€ B w. (27) 

Here w includes the energy of baryons, radiation, and e^; it also includes the neutrino energy if the 
disk is i/-opaque. The field has the strongest effect on light charged particles — e^. It introduces 
the discrete energy levels (Landau & Lifshitz 1980), 

1/2 

J = 0,1,.., (28) 



l+pl + 2 3 — I 



m e <r 

where — oo < p z < oo is the component of the electron momentum parallel to the field and 
Ub = eB/m e c. The magnetic field also changes the phase-space factor in equation (1): p 2 dp 
is replaced by [MIb / 'm e c 2 )(dp z /2) . Both effects are important if {ft£lB/m e (?) > p 2 z - The mean 
parallel momentum of the relativistic e ± equals y/3kT/c, and the condition for the field effects to 
be important reads 

> 3^ 2 . (29) 

For any plausible e# < 1, the magnetic field is important only where the electrons are degenerate, 
and the energy density w is dominated by either baryons, Wb = (3/2)kTp/m p , or degenerate 
electrons, w^eg = (9/4)(7r/3) 4 / 3 /zc(y e /9/m p ) 4 / 3 . The condition (29) can be written as 

d< \m p ) \ (vr/3) 7 /i2 (6a/es) i/4 ye i/3 €B K (7r 3 F 4 /32a/) _ m 
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where ctf is the fine structure constant. The upper line corresponds to Wf, > Wd cg , and the lower 
line to Wdeg > Wb- 



2.2. GRB central engines 

There exists a general constraint on the electron degeneracy in the engine. Its derivation 
makes use of two facts: (1) The engine is gravitationally bound — otherwise it would explode, and 
a high baryon contamination of the fireball would be inevitable. GRB models normally envision a 
quasi-static engine that liberates a fraction of its gravitational binding energy and passes it to a 
tiny amount of mass outside the engine, thus creating a highly relativistic outflow. (2) The engine 
is compact (size r < 10 7 cm), and it has a relativistic blackbody temperature 9 = kT/m e c 2 > 1. 

The sound speed in a gravitationally bound object of size r and mass M must be smaller than 
(GM/r) 1 / 2 . It gives the constraint 

7 S0 ' lc2 (£)"'• <31> 

where r g = 2GM/c 2 is the gravitational radius of the object. Pressure P = P y + P± + P v + Pj, 
includes contributions from radiation, e ± , neutrinos, and baryons. We have P 1 + P± = (ll/12)aT 4 
if the are weakly degenerate, and a maximum P v = (7/24) aT 4 for each neutrino species if it 
is thermalized, where a = (TT 2 k 4 /15h 3 c 3 ) = 7.56 x 10~ 15 erg cm~ 3 K~ 4 is the radiation constant. 
Approximately, 

P w aT 4 + -P-hT. (32) 

When comparing the two terms in equation (32), it easy to see that 

aT 4 7T 2 9 3 



(p/m p )kT 15 03 



(33) 



dcg 



Hence, the baryonic pressure gets dominant at 9 ~ #deg [and at even lower temperatures, 9 < 
(Tr/4)Ye^ 3 9f\ es , the pressure of degenerate electrons Pdcg = (3/4)(-7r/3) 4 / 3 7ic(Y e p/mp) 4 / 3 takes over]. 
Our goal is to derive an upper bound on #/6*d C g and therefore we consider 9 ^ #deg with P k, aT 4 . 
Constraint (31) then reads 

which also gives a lower bound on the electron chemical potential (eq. 13), 

We know from § 2.1 that the condition Y e < 0.5 reads fi > Q/2 for ^-transparent matter and fi > Q 
for z/-opaque matter. We now see that any engine of size r and density p > 10 10 (r/3r 9 ) _1 g/cm 3 
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satisfies this condition and tends to an equilibrium Y e < 0.5. This is evidently the case in models of 
neutron star mergers (Ruffert & Janka 1999), as well as magnetized neutron stars. The collapsar 
scenario (MacFadyen & Woosley 1999) invokes a relatively low-density accretion flow, and here 
neutronization is questionable. We therefore shall study accretion flows in more detail. We also 
need to check whether the equilibrium Y e is achieved on the accretion timescale. 



2.2.1. GRB accretion flows 

All accretion models of GRBs invoke rotation that creates a funnel along which the fireball 
can escape. The accretion flow is viewed as a rotating disk maintained in hydrostatic balance in 
the vertical direction, which gradually spirals to the central black hole. A standard model assumes 
turbulent viscosity in the disk with a stress tensor W r< ^ = ac 2 , where c s = (Pj p) 1 ! 2 is the isothermal 
sound speed and a = 0.01 — 0.1 (Balbus & Hawley 1998). Same stress tensor can be described in 
terms of a viscosity coefficient v = (2/3)ac s H, where H is the half-thickness of the disk. 

The standard disk theory gives the velocity of accretion at radius r 

W r<t> or \ i fiiA 1 / 2 /q£\ 

u = r^s' s{r) = 1 -{V) ' (36) 

where J)k(?") = (GM/r 3 ) 1 ' 2 is the angular velocity of Keplerian rotation, M is the black hole mass, 
and ri n is the inner boundary of the disk (the marginally stable orbit); r- m = 3r g for non-rotating 
and r; n = r g /2 for extremely rotating black holes. 6 S(r) varies from to 1, and it equals 0.5 at 
a characteristic radius r = 4r; n where the integrated dissipation rate peaks (AirHr 2 q + oc S/r, see 
eq. [40] below). One expects rotating black holes in GRBs, and r rs 3r g is a reasonable characteristic 
radius. 

The vertical hydrostatic balance reads c s = HQk, which we use to write the accretion time as 

, , r S fH\- 2 - f2H\' 2 f a \~W r \ 3/2 / M \ , s 

( « W = ^ = ^(7j =2.9X10-S( V ) (_) (-) (37) 

We neglect here the disk gravity, which is a valid approximation as long as Mt a <C M where M is 
the accretion rate. The typical t a is much shorter than the burst duration, and accretion is viewed 
as a quasi-steady process. It can power a relativistic outflow ("fireball") with luminosity 



where e/ is the efficiency of Mc 2 conversion into a fireball, which is below the net efficiency of 
accretion (e ~ 0.1). Most of the accretion energy is either carried away by neutrinos or advected 
by the accretion flow (in case of small neutrino losses). 



6 We give here simple Newtonian estimates and trace the hole spin only through its effect on n n ; other relativistic 
corrections are modest and weakly affect the results. 
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The disk baryonic density is given by 

Mt a nn ininn • (2H\^ ( r \~ 3/2 ( a \ -i / M 



47rr 2 iJ 

and it is heated viscously with rate 

. 3M^ f2H\~ 1 ( r \ -4 / M x ~ 3 



8ttH 



^ ps 5.1 x 10 33 SM 32 (J^j ergcm-V 1 . (40) 



The flow has a huge optical depth for radiation and the radiation is trapped — its diffusion is 
negligible on the accretion timescale. The only cooling mechanism of the flow is neutrino emission 
which becomes efficient at M ^ 10 31 g s _1 and small r (e.g. Popham, Woosley, & Fryer 1999; 
Narayan, Piran, & Kumar 2001; Kohri & Mineshige 2002). An upper bound on the temperature 
is derived from the assumption that the neutrino cooling is absent. Then the flow does not loose 
the dissipated energy and instead traps it and advects. Its energy density can be estimated as 
io miH ps q + t a ps (3/8)(r g /r)pc 2 (we use eqs. [37] and [40] with S = 0.5). The internal energy of such 
a hot advective flow is dominated by radiation and e^, so 

11 3r 

«W = jar max ps -gV 2 , (41) 

where factor 11/4 accounts for the contribution of relativistic weakly degenerate (neutrino 
contribution can further increase this factor if the neutrinos are reabsorbed, and then T max will be 
slightly lower). Equation (41) yields 

/ \ -1/4 

1/4 / r • 



kT max ps 13^ l^—j MeV. (42) 

The actual temperature can be significantly lower if the neutrino cooling is significant. The main 
cooling process is capture on nucleons: e~ + p — > n + v and e + + n — > p + i>, which also shapes 
Y e as we discussed in § 2.1. We now evaluate the characteristic accretion rate M eq above which the 
e ± capture is rapid enough to establish an equilibrium Y e . 

The equilibrium Y e is achieved when the flow has emitted one neutrino per nucleon. It is easy 
to see that disks with efficient neutrino cooling always reach the equilibrium. Indeed, the mean 
energy of the emitted neutrinos, E v ^ 5kT, is below the liberated accretion energy per nucleon, 
E n = 100 — 300 MeV, and an efficient cooling implies that more than one neutrino per nucleon 
is produced. Thus, M cq should be looked for in the inefficient (advective) regime with T ps T max . 
Such a flow is only mildly degenerate, and the rates of e ± capture read (in zero order in fi/9) 

ii e - p ps 1.5 x 10~ 2 n p 6 5 cm~V\ (43) 
h e+n ps 1.5 x 10- 2 n n 6 5 cnrV 1 . (44) 

The neutronization timescale is t n = n p /h e - p ps 70#~ 5 s, and it should be compared with the 
accretion timescale t a (eq. 37), 

^™^(-V- («) 
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We substitute T = T max , take into account the hydrostatic balance 

H ( w max \ 1/2 1 1 



r \ 3p / tl^r 2' 
and use equation (39) with S = 0.5 to get 



(46) 



We conclude that disks with 



o n / r \ 13/10 / a \ 9/5 / M \ 6/5 
M>M„ = 3.8xlO-(-) (-) (-) gs- (48) 

have t n < t a and hence approach the equilibrium Y e . Radius r entered this expression in the 13/10 
power, which implies that the characteristic M eq depends significantly on the black hole spin: the 
characteristic radius (where viscous dissipation peaks) decreases from r ~ 10r s to r ~ r g as the 
hole spin increases from zero to a maximum value. 

The equilibrium Y e is below 0.5 if the flow temperature is below T n that was calculated in § 2.1 
(eqs. [14] and [23]). It is instructive to compare T n with the maximum accretion temperature T max 
(eq. 42), 



± max \ «' g 



where k = 2.5 for ^-transparent and k = 1.8 for z/-opaque flows. The condition for a neutron- 
excess, T < T n , is most difficult to satisfy in low-M flows where the neutrino cooling is inefficient 
and T w T max . Such flows are zv-transparent, so that n = 2.5 should be used in equation (49). 
Substituting equation (39), T = T max , and H/r = 1/2, we find that T max < T n if 

M>M^7.6x 1 030 (^) 1/2 (^)(^-) 2 g8 -. (50) 

The "neutronization" accretion rate M n is comparable to M cq . Plausible M in GRB accretion flows 
are 10 32 g/s and higher, and they should have a neutron excess. 



2.3. De-neutronization in the fireball? 

Are outflows from neutron-rich engines also neutron-rich? The fireball picks up baryons from 
the surface of the central engine, and the surface density is relatively low. Y e might change while 
the matter is escaping into a fireball. 

Let us consider fireballs produced by accretion disks. The disk is turbulent (the turbulence is 
the source of viscosity that is responsible for accretion), and it is mixed in the vertical direction 



-14- 



on the sound -crossing timescale t m i x « H/c s (the turbulent velocity is somewhat smaller than c s , 
however, the thickness of surface layers is also smaller than H, and t m i x ~ H/c s is about right). 
Given the hydrostatic balance, H = c s /Qk, one gets 

The turbulent material circulates rapidly up to the surface and back to the interior of the disk, 
and a small portion of it can (also rapidly) escape in each circulation. The escape timescale of the 
fireball is comparable to Sl^ 1 . As an element of matter elevates to the surface its Y e would increase 
if it adjusted instantaneously to a new equilibrium value. However, the time fi^ 1 the element has 
before it sinks back into the disk (or escapes) can be too short for the adjustment. Then Y e of 
the escaping matter corresponds to p and T inside the disk, where it has spent almost all the time 
before the sudden escape. 

To check this picture, let us evaluate the timescale of "de-neutronization" of an initially neutron 
rich material that has suddenly expanded into a low-density, hot fireball. Neutrons tend to convert 
back into protons via two charged current reactions: e + capture and v absorption (/3-decay is 
slow and negligible). The fireball has temperature 9 > 1 and non-degenerate electrons, \x <C 9. 
Equation (44) then gives the e + capture timescale 

' nn ~ 70 B, (52) 



n. 



<:~ n 



9 5 



For fireball temperatures up to 8 MeV the e + -capture is slow compared to f2 K x and does not affect 
Y e of the escaping material. 

The v absorption by neutrons is potentially more important for de-neutronization. The cross 
section of this reaction is 

1/2 



a a (u) = 2.4 x 1(T 44 [1 -f-(u + Q)] (u + Q) 2 , . { + Q) , 
and the corresponding rate of de-neutronization is 



1 

1 



cm 2 w 2.4 x 1(T 4 V cm^, (54) 



h n = -cn n I n^cjaduj, (55) 
J o 

where = dn u /du> is the spectrum of the neutrino number density. This gives the timescale of 
de-neutronization, 

n n 1.4 x 10 33 

tv = tt-t = ==— s, (56) 

\n n \ n u uj z 

where bar signifies an average over a neutrino spectrum. An upper bound on the neutrino density 
outside the disk is given by 

n « < 8^' (57) 
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which states that the energy flux of the (electron) neutrino from the two faces of the disk, 2F U m 
2n u uJm e c 3 / 7r , is smaller than the accretion energy released per unit time per unit area of the disk, 
F = (3/47t)GMMS/r 3 . Thus, we get (with S ~ 0.5) 

3/2 



^K>10 2 =M 3 - 2 M— -) . (58) 



^ 32 y M& ) y 3rg 

The mean uj and uj 2 are determined by the state of material that emits the neutrinos. If the disk 
is zv-transparent, the neutrinos have the spectrum (see eq. [6]) 

oc (uj + QY 



uj + Q- H 
exp ( - Q ) + 1 



1 



nl/2 






2 ^ 4 
UJ PS UJ 


exp (- 





1 , (59) 



where is temperature inside the disk, and Q and ^/ have been neglected compared to the typical 
uj = uj. This gives uj = (31/6)[C(6)/C(5)]0 » 5.070 and = (63/2)[C(7)/C(5)]# 2 « 3O.60 2 . 

The optical depth of the disk for neutrino absorption by neutrons is 7 

Here we assumed T > Td cg , so that Y e is not much below 0.5 and 2H/r = (T/T max ) 2 (P oc T 4 , 
sec cqs. [32] and [33]). From equations (58) and (60) one can see that ^-transparent disks (r„ < 1) 
have > 1 for a > 0.01, and hence their fireballs do not experience any significant absorption 

of neutrinos. 

High-M disks with neutrino luminosity L v > 10 53 erg/s create a sufficiently dense bath of 
neutrinos that can impact Y e of the ejected fireball. High-M disks are z/-opaque and the emitted 
neutrinos can be approximated by a blackbody spectrum with uj ~ (7/2)[^(4)/C(3)]0^ ~ 3.15#j, 
and uj 2 « 15[C(5)/C(3)]03 ~ 12. 96* 2 , where 6 V is the temperature of the ^-photosphere in the disk. 
The disk also emits v from a corresponding ^-photosphere. Both v absorption by neutrons and v 
absorption by protons occur at the base of the fireball and set a new equilibrium Y e such that the 
rates of v and v absorptions are equal. One can expect the new Y e to be below 0.5 (see Qian et al. 
1993, where a similar problem is discussed in the context of supernova engines, and Pruet, Fuller, 
& Cardall 2001). This expectation is based on two facts: (1) the number densities of v and v are 
approximately equal, \n v — n D \ <C n„, — otherwise they carry away too large a leptonic number 
from the disk, and (2) the absorption cross section is proportional to neutrino energy squared (in 
the limit uj j> Q, which is valid for the hot high-M disks). It gives the equilibrium n n /n p w Op/0 2 . 
where V and Op are temperatures of the v and v photospheres. One expects V < Op because the 
i^-photosphere is likely closer to the disk surface (the neutronized disk is more opaque for u) and 
then n n /n p > 1. This argument assumes, however, that the disk temperature decreases toward its 
surface. It may not hold if the dissipation rate peaks at the surface and makes it hotter than the 
interior of the disk. 



Absorption makes a major contribution to the disk opacity for neutrinos. 
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3. Nucleosynthesis 

The ejected fireball has a low baryon loading and a high temperature, and its nucleons are 
initially in the free state. As expansion proceeds, the fireball cools adiabatically and, when its 
temperature decreases to kT* ~ 100 keV, fusion reactions shape the nuclear composition like they 
do during the primordial nucleosynthesis in the Universe (Wagoner, Fowler, & Hoyle 1967). In both 
cases, we deal with an expanding blackbody fireball with initially free nucleons, and the following 
three parameters control the outcome of nucleosynthesis: (1) photon-to-baryon ratio (f> = n 7 /nft, 
(2) expansion timescale r* at the time of nucleosynthesis, and (3) n/p-ratio prior to the onset of 
nucleosynthesis. In the Universe, 4> « 3 x 10 9 , t* k, 10 2 s, and n n /n p « 1/7. Below we formulate 
the nucleosynthesis problem in GRBs and give a qualitative comparison with the big bang. Then 
we make detailed calculations of nuclear reactions in GRB fireballs. 

3.1. Fireball model 

A customary GRB model envisions a central engine that ejects baryonic matter at a rate 
Mb [g/ s ]j thermal energy at a rate L t h 3> M^c 2 , and magnetic energy at a rate Lp. The produced 
fireball expands and accelerates as its internal energy is converted into bulk kinetic energy The 
fireball is likely to carry strong magnetic fields and the Poynting luminosity Lp may be higher than 
the thermalized luminosity. In this case, the magnetic fields prolong the acceleration after all the 
thermal energy has been converted into bulk expansion. A maximum Lorentz factor of the fireball 
can be estimated as r max = (L t h + Lp)/M&c 2 , and it is at least 10 2 for GRB explosions. 

The nucleosynthesis occurs when the fireball temperature drops to about fcT* ~ 100 keV at 
a radius R* ~ 10 7 — 10 9 cm. The timescale of expansion to this radius, R*/c, is shorter than the 
duration of the engine activity. Therefore, in the nucleosynthesis calculations, the fireball can be 
modeled as a quasi-steady outflow. 

Let the outflow expand in an axisymmetric funnel with a cross section 

S(R) = S (^y. (61) 

For example, ip = 2 for a radial funnel (and also for a spherically symmetric explosion), and ip = 1 
for a parabolic funnel which may develop in a collapsing progenitor of the GRB (Meszaros & Rees 
2001). The outflow is a relativistic ideal fluid with baryon density p, pressure P, and energy density 
w = 3P 3> pc 2 ; all these magnitudes are measured in the fluid rest frame. In spherical coordinates 
x l = (t,R,6,<p) the outflow has 4-velocity u % = dx l /dr = (u* ,u R ,u ,u^), where r is proper time. 
We assume = and Ru e <C u l ,u R . The latter assumption is satisfied at all radii for a radial 
explosion (u e = 0) and at R » vq for a collimated explosion (ip < 2). 

The outflow dynamics is governed by the conservation laws Vi(pu l ) = and Vj(T^) = 0, where 
T % k = u % Uk(w + P) + P5 l k is the stress-energy tensor. The electromagnetic tensor is not included here, 
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which greatly simplifies the problem (we are interested in the early hot stage when the expansion 
is likely driven by the thermal pressure even in the presence of strong fields). Then the baryon and 
energy conservation laws read 

Spu R = M b , S{w + P)u l u R = L th . (62) 

The high Lorentz factor of the expansion r = > 1 implies u R /c « T. Equations (62) then yield 

Mb w 3 Lt h 



STc' p 4M b T' 



(63) 



We assume that the outflow does not exchange mass or energy with the surroundings, i.e., M b (R) = 
const and L t h(R) = const. Equation (63) then gives w oc p/T while the first law of thermodynamics 
d(w/p) = —Pd(l/p) gives w oc p 4 / 3 . Excluding p from these two relations one gets w 1 ^ oc 1/T and 

TV = T , (64) 

where To is a constant. Equation (64) is strictly valid at temperatures kT <C m e c 2 = 511 keV 
where the e ± energy density can be neglected. At small T (where kT J> m e c 2 ), the inclusion of 
reduces T by a modest factor (11/4) -1 / 4 . 

Equation (62) implies w = aT 4 oc (r 2 ^) -1 and, combining with equation (64), one gets 

The outflow may be collimated already at its base into a solid angle f^o = SbAo < ^7r — its 
transonic dynamics near the central engine is complicated and unknown. With reasonable accuracy, 
the simple estimate L t h ~ SqclTqC gives 

The value of T is most sensitive to the engine size r . It can be as small as 3 x 10 5 cm (if the 
outflow is powered by a Kerr black hole via the Blandford-Znajek process) or as large as 10 7 cm (if 
the outflow is powered by an accretion disk). 



r-%-(°y -(*)*". (65) 



A major parameter of the nucleosynthesis problem is the ratio of photon density, 

_ 30C(3)u> w 



(67) 



7 7T 4 fcT 2.70&T' 

to baryon density n b . This ratio is constant during the expansion, 8 and its value is (cf. eqs. [63] 
and [64]) 

j. _ n 7 _ 3 m P ^th - - - ~/i / nth \ ( kTo \ ~ ' L t 



n b 4 M b 2.7kT 



7 7 „ 1 n4 ( Vth \ ( kT \ 1 L th 



8 It is proportional to entropy per baryon: s/k = (27r 4 /45C[3])n 7 /n(, w 3.602n 7 /n(, where k is Boltzmann constant, 
and it must be constant in an adiabatic expansion. 
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The typical is 5 orders of magnitude smaller compared to that of the Universe. 



The nucleosynthesis must occur during the acceleration stage of the fireball because the accel- 
eration ends at a low temperature T = T /r max <C T*. The accelerated expansion is described in 
the comoving time by dr = dt/T ps dR/cT, which gives (we use eq. [65]) 

■log©) 



r(T) = 



ro 
c 



la. 



r x(V-2)/V 



Y> = 2, 
V> < 2. 



(69) 



The time of nucleosynthesis depends on the shape of the funnel: r(i2*) ~ 3ro/c for ip = 2 and 
t(R*) ps 2(To/r*)(ro/c) for ip = 1. The timescale of the density fall off is 3 times shorter than that 
of the temperature because p oc T 3 . Therefore a reasonable choice of the characteristic expansion 
timescale during nucleosynthesis is t* ps ro /c for ^ = 2 and t* ps (2ro/3c)(To/T,,,) for ip = 1. 

Note that L t h, and To enter the nucleosynthesis problem only in combinations that deter- 
mine <fi and t(T), and play no other role. For example, the initial temperature is not important as 
long as kT > 200 keV. 



3.2. Simple estimates and comparison to the big bang 

First we evaluate the temperature T* at which we expect the nucleosynthesis to begin. Let 
us remind why kT* ps 80 keV in the big bang. If the nuclear statistical equilibrium (NSE) were 
maintained, recombination of nucleons into a-particles would occur at kT ~ 200 keV (e.g. Meyer 
1994). The reason of the nucleosynthesis delay until kT ~ 80 keV is what is sometimes called 
"deuterium bottleneck". Before fusing into a-particles, the nucleons have to form lighter nuclei 
- deuterium, tritium, or 3 He. At kT > 100 keV, these elements have very low equilibrium 
abundances, which implies a very long timescale for their fusion, and helium is not formed even 
though it is favored by the NSE. In particular, the deuterium abundance is suppressed by the very 
fast photodisintegration 7 + d — > n + p that balances the opposite reaction n + p^d + ^dX the 
equilibrium value 

Y d « 7 x lO-^YnYpT^ 2 exp , (70) 

(e.g. Esmailzadeh, Starkman, & Dimopoulos 1991). Here Yi = rii/rib is abundance of the i-th 
element and Tg = T/10 9 K = kT/86.17 keV. Nucleosynthesis starts when the exponential wins the 
pre-exponential factor to give a noticeable Yd ~ 10 -3 . Substituting cj) = 3 x 10 9 , one finds that it 
happens at kT* ~ 80 keV. 

In GRB outflows, is 5 orders of magnitude smaller, and the NSE favors the recombination 
of nucleons into a-particles at a temperature as high as 500 keV. For the same reason as in the 
big bang, the nucleosynthesis is delayed to a lower T*. Using equation (70) with (f> = 10 5 , we find 
kT* ps 140 keV, which is almost twice as high as the big bang T* — it depends logarithmically on 
0. 
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At a given T, nuclear reaction rates scale as Y oc p oc (here Y is an element abundance). 
The ratio of the expansion timescale r* to a reaction timescale, T rea c = Y/Y, behaves as T*/<f). This 
combination is ~ 3 — 30 times smaller in GRBs compared to the big bang, which may seem not a 
crucial difference. In fact, the big bang nucleosynthesis was dangerously close to a freezeout of the 
neutron capture reaction n+p — ► d + 7, without which nucleosynthesis cannot start. Therefore the 
difference by a factor of 10 can be crucial and it is instructive to compare accurately the capture 
timescale with the expansion timescale in GRBs. 

The neutron capture rate varies slowly with temperature (see Fig. 11 in Smith, Kawano, Sz 
Malaney 1993), and near 100 keV it is approximately Y c ps 2.5 x 10 4 pY p Y n s _1 . The reaction 
timescale is 

min{y p ,y n } _ 1.2 x 10-y 

Y c T| max{F n , Y p } S ' ^ 

It should be compared with the timescale of the density fall off, r*. For the big bang, r*/r c « 30. 
For GRBs, we get 

re « 10- 4 (4 V£ V 1 s, (72) 



10 5 / V 0.5 



— W l) V> = 2, 



Z* ~ J V 105 / V 3xl ° 6 / V0 - u/ (7^ 



W I \ 3xlO e 



One can see that, in a radial explosion (-0 = 2), the neutron capture rate is marginal for a successful 
nucleosynthesis. In a collimated explosion, the ratio t*/t c is higher by a factor of To/T* ~ 10, and 
nucleosynthesis is efficient. 

An important difference between the GRBs and the big bang is the n/p-ratio. In the big bang, 
n n /n p = 1/7 (Y e = 7/8), which leads to 25% mass fraction of helium after the n-p recombination, 
while 75% of mass remains in protons (and a tiny amount of other nuclei). In GRBs, n n /n p >1 
(Y e < 0.5), and there are leftover neutrons even if all protons are consumed by helium production. 
The minimum mass fraction of leftover neutrons is 

X n = 1 - 2Y e . (74) 



3.3. Detailed calculation 

3.3.1. The code 

We will keep track of elements with mass numbers less than 5 (like the big bang, the abundances 
of heavier nuclei are very small). The six elements under consideration are neutrons, protons, 2 H, 
3 H, 3 He, and 4 He; they are denoted by n,p,d,t,3,a, respectively, and the photons are denoted by 
7. Their abundances are measured by Yj, = ni/n^, or mass fraction = AjXi, where n« and A{ are 
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the number density and the mass number of the i-th species. The photon and matter densities are 
known functions of temperature, 

2 fkT\ 3 

n 7 = -2 C(3) ( — J = 2.02 x 10 28 T| cm" 3 , p = 3.39 x W 4 ^' 1 ^ g cm" 3 . (75) 

The evolution of nuclear composition is described by the set of equations 

Yi = ^YkYAklifi-^YiYjlijkl], (76) 

Y n = -Y d -2Y t -Y 3 -2Y a -^, (77) 

Tp 

Y p = -Y d -Y t -2Y 3 -2Y a + ^. (78) 

Here dot signifies a derivative with respect to proper time r, [ijkl] = UbWvij^ki is the rate of 
reaction i + j — > k + I, and all quantities are measured in the rest frame of the outflow. It is 
sufficient to calculate the reaction rates for deuterium, tritium, and helium isotopes 3 He and 4 He 
(i = d, t, 3, a in eq. 76). Then Y n and Y p are found from the neutron and proton conservation laws 
(eqs. [77] and [78]). The latter include the Yn/rp term — the conversion of neutrons into protons 
via /3-decay with the mean lifetime of neutrons rg = 900 s. The /3-decay is negligible in GRBs and 
we keep it for the code tests on big bang nucleosynthesis. 

The sums in equation (76) are taken over all possible reactions with participation of the i-th 
nuclei. Not all reactions are important (Smith et al. 1993). For example, reactions that destroy 
a-particles can be neglected as Y a is far below its equilibrium value. We include the following 
reactions in the calculations: n + p <-> d + 7, n + 3 <-> p + t, n + 3 — > a + 7, p + d ^ 3 + 7, 
p + t — > a + 7, d + d — >■ p + t, d + d — > n + 3, (i + t^n + a, d + 3^p + a, and take their rates 
from Smith et al. (1993) and Esmailzadeh et al. (1991). We then have 

Y d = Y n Y p [npd 7 ] - Y d Y 7 [d 7 np] - Y p Y d [pd3 7 ] - 2Y d 2 [ddn3] - Y d Y 3 [d3pa} 

-Y d Y t [dtna] - 2Yj [ddpt] , (79) 

Y t = Yj [ddpt] + Y n Y 3 [n3pt] - Y d Y t [dtna] - Y p Y t [ptn3] - Y p Y t [pto 7 ] , (80) 

^3 = Y p Y d [pd3 7 ] + Yd 2 [ddn3] + Y p Y t [ptn3] - Y d Y 3 [d3pa] - Y n Y 3 ([n3pt] + [n3a 7 ]), (81) 

Y a = Y d Y t [dtna}+Y d Y 3 [d3pa}+Y n Y 3 [n3a~f}+Y p Y t [pta 7 }. (82) 

The reaction rates are functions of T and p, and the set of equations is closed by equation (69) that 
relates T and r. We solve equations (77-82) numerically with initial conditions Y p = Y e , Y® = 1—Y e , 
Y a = 0, and Y d = Yt = Y 3 = at an initial temperature kTo ~ 1 MeV. At high temperatures, 
the abundances of all elements except n and p are negligibly small. The abundances of deuterium, 
tritium, and 3 He are close to quasi-steady equilibrium at kT > 150 keV, i.e., the rates of their 
production Y + and sink Y~ are much higher than the expansion rate Y/t, and hence Y + and Y~ 
almost balance each other. Where an element abundance Y approaches the equilibrium value, it is 
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calculated by setting the net Y = Y + — Y~ = 0. (Thus, one avoids numerical integration where Y 
is a small difference of big numbers.) The code has been tested with the big bang nucleosynthesis. 
It successfully reproduced the standard evolution of the cosmological nuclear composition. 



3.3.2. Results 

Figure 3 compares the nucleosynthesis in GRBs and the big bang (BB). Parameters of the 
GRB fireball in this example are chosen as (j> = 10 5 , vq = 3 x 10 6 cm, ip = 2, and Y e = 0.5. 
As expected, GRB nucleosynthesis occurs at higher temperatures. The freezeout mass fractions 
of helium and deuterium are X a rs 0.16 and « 0.03, and about 81% of mass remains in free 
nucleons. Interestingly, the deuterium evolution is qualitatively different in GRBs: X& increases 
monotonically and freezes out at the 3% level which is almost 3 orders of magnitude higher than 
in the big bang (see also Pruet et al. 2002). This difference is due to a high abundance of neutrons 
and a negligible /3-decay. As a result, the deuterium production by the neutron capture reaction 
outweighes its burning rate and X& grows. Also, the freezeout happens quickly in the radial 
expansion (T oc exp[r/ro] where tq = r$/c) and Xd could not decrease much below 1% even if the 
neutron capture reaction were completely switched off at, e.g., kT < 80 keV. 

The post-nucleosynthesis fireball is predominantly composed of free nucleons and helium, and 
it would be useful to know how their freezeout ratio 

Ji. n -+■ A. p 

depends on the four parameters of the explosion ro, <f>, Y e , and ip. Figure 4 shows f(ro,<f>) for 
radial explosions with Y e = 0.5. Contours / = const on the ro — <p plane are perfect straight lines 
<p oc Tq -0 ' 85 (or ro oc (f)^ 1 ' 18 ). This implies that / depends on combination ro^ -1 ' 18 . A similar 
situation takes place for collimated explosions (Fig. 5). The slope of all contours in Figure 5 is 0.90 
and / can be viewed as a function of ro0~ L11 . Calculations with Y e / 0.5 give different values of 
/, however, contours /(ro, <p) = const have exactly the same slope as in the case of Y e = 0.5. 

The results of nucleosynthesis calculations can be represented in a compact form if we define 
a new variable 



Figure 6 shows the outcome of nucleosynthesis /(£) for different Y e . Strikingly, /(£) is practically 
identical for radial and collimated explosions. Figure 6 gives a full description of helium production 
in relativistic explosions: from this figure and equation (84) one finds / for given 4>, ro, Y e , and 
Y> = 1,2. 

One can see that nucleosynthesis is most efficient in Y e = 0.5 models, where X® = In 
the most interesting region of parameters, where / changes from / < 1 to / > 1, one observes a 
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Fig. 3. — Evolution of deuterium and helium abundances with temperature in an expanding fireball. 
The GRB case is shown by thicker curves for a radial explosion with Y e = 0.5, (ft = n ^/ n b = 10 5 , 
and assuming a central engine of size ro = 3 x 10 6 cm (tq sets the expansion timescale, cf. eq. [69]). 
For comparison, the big bang (BB) nucleosynthesis is also shown (with (ft = 3 x 10 9 ). The dotted 
curves display the equilibrium (production=sink) abundances of deuterium. 
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Fig. 4. — Contours of freezeout ratio / = X a /(X p -\-X n ) on the rQ—<p plane for radial explosions with 
Y e = 0.5. Right axis shows (L^/ Mt ) c 2 )(kTo/MeV)~ 1 , which is related to </> by equation (68). If the 
fireball is not Poynting-flux dominated, Lp < L t h, its final Lorentz factor equals r^h = L^/M^c 2 . 



-24- 




Fig. 5. — Same as in Figure 4 but for parabolically collimated explosions. 
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simple behavior 

/ = £, 0.03 < / < 30, (85) 

for both radial and collimated explosions (this explains our choice of constants 0.2 and 7.6 in 
eq. [84]). 

Fireballs with Y e / 0.5 have an upper bound on the freezeout abundance of helium maxl a = 
2min{y e , 1 — Y e }. For explosions with a neutron excess, Y e < 0.5, this gives a maximum /, 

2Y 

/max = i_2Y ' 

With increasing £, the dependence / « £ switches to / = / max = const (Fig. 6). 

The post-nucleosynthesis fireball is dominated by a-particles if / > 1, which requires Y e > 0.25 
and £ > 1. For example, in a radial explosion with Y e = 0.5 it requires 4> < 2.5 x 10 4 (ro/3 x 10 6 ) - 85 , 
and a similar condition for parabolic explosions reads <f) < 6.2 x 10 4 (ro/3 x 10 6 ) ' 9 . The typical 
GRB parameters happen to be just marginal for a successful nucleosynthesis: / varies from 10~ 2 
to / m ax i n the expected range of 4> and r . 



4. Spallation 

Synthesized helium may be destroyed during the subsequent evolution of the explosion. Spal- 
lation reactions can occur when an a-particle collides with another particle with a relative energy 
exceeding the nuclear binding energy. The fireball temperature at the acceleration stage is too low 
for such reactions, however, the collisions may become energetic if (1) there appears a substantial 
relative bulk velocity between the neutron and ion components or (2) internal shocks occur and 
heat the ion fireball to a high temperature, much above the blackbody value. 



4.1. Neutron-ion collisions during the acceleration stage 

Near the central engine, the neutron component of the fireball is well coupled to the ion 
component by elastic collisions with a small relative velocity (3 = v/c (Derishev et al. 1999, 
Bahcall & Meszaros 2000). The collision cross section can be approximated as U{ = where 
uo ~ 3 x 10 -26 cm 2 if the ions are protons (and oo approaches 10~ 25 cm 2 for a-particles). In the 
fluid frame, the mean collisional time for neutrons is 

11 

Tcoii = = , (87) 



where n« is the ion density. The ion density is a fraction of the total nucleon density rib which 
behaves as = const /T 3 during the acceleration stage (see eqs. [63] and [65]). The constant of 
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Fig. 6. — Outcome of nucleosynthesis in relativistic explosions for Y e = 0.1, 0.25, 0.4, and 0.5. The 
ratio of mass fractions of helium and free nucleons, / = X a /(X n + X p ), is shown as a function of 
£ defined in equation (84). Dotted and solid curves display radial (ip = 2) and collimated (ip = 1) 
explosions, respectively. Dashed line shows / = £, which well approximates the numerical results 
in the range 0.03 < f < /max (if Ye is not much below 0.5). 
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proportionality can be expressed in terms of r/ t h = L^/M^c and 4> = n^/rib, 

n b = m 3 A = 5.61 x 10 38 ^ cm -3. (88) 

4vr 14 V ft 7 </> 4 r 3 </> 4 r 3 v ; 

In the fixed lab frame, the ion fluid is accelerated by radiative or magnetic pressure, and its Lorentz 
factor is doubled on timescale t = R/c. Neutrons "miss" the ion acceleration by AT/F ps i co iiA = 
Tct co \\/R < 1 and have a smaller Lorentz factor r n = T — AT. The relative velocity of the neutron 
and ion components is (3 = (T 2 - Tl)/(T 2 + T 2 n ) « (r - T n )/F, which is 

t rLniCjQ R 

The energy of neutron-ion collisions becomes sufficient for spallation reactions if (3 exceeds 
f3 sp ~ 0.1. This can happen at late stages of the fireball acceleration, at high T but not exceeding 
Tmax = L/MfrC 2 = r?th + VP (L = L t h + Lp is the total luminosity of the fireball that includes the 
Poynting flux). In the case of T(R) « R/tq (radial explosion), (3 reaches (3 sp when V reaches 

r -^(¥)(^)" 4/3 fe) 1/3 fe) 1/3 (S 1/3 ^ 



Spallation takes place if r sp < r max , which requires 4> > 4>sp, 

3/4 



<Asp « 4 X 10 4 (^) . 



(91) 



We keep here L t h/L < 1 as the most uncertain parameter which may be much below unity. Note 
also, that in the case of a non-radial explosion, ro should be replaced by R/Y in equation (90). 
Using equation (68) it is easy to show that the condition <p > 4> sp is equivalent to 

G^x^r <92> 

Let Y n be abundance of neutrons that survived the nucleosynthesis, and suppose that <j> > 4> sp . 
The lifetime of a-particles bombarded by neutrons with (3 ps (3 sp is iiif e w {t co \\/Y n ){<jQ / o SVl ) 
(Psp/Yn)t, where o" sp ~ o"o is the spallation cross section (it is roughly equal to the geometrical size 
of the nucleus -nr 2 N 5.3 x 10~ 26 A 2 / 3 cm 2 with A = 4 for a-particles). Hence a modest Y n ~ 0.1 
should be sufficient for a significant spallation. 

In case of Sj V6ry low bciryon loading (high r max ) the neutrons decouple from the ions before 
the end of the acceleration stage and their Lorentz factor T n saturates at T^cc < r max . The relative 
velocity (3 approaches unity when the decoupling happens, and the last neutron-ion collisions arc 
energetic enough for pion production which leads to multi-GeV neutrino emission (see Derishev 
et al. 1999, Bahcall & Meszaros 2000). The fireball Lorentz factor at this moment, Tdec, is given 
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by equation (90) with (3 sp replaced by unity, i.e. r^c = 10 1//3 r sp . The decoupling takes place if 
r<iec < r max , which requires 4> greater than 

0dec«7xlO 4 (^ 3/4 . (93) 

The condition <f> > </>dec is equivalent to r max > 300(/cTo/MeV)(L/L t h) 1//4 . The decoupling is always 
preceded by spallation. 

The upper bound on the neutron Lorentz factor due to decoupling is essentially determined by 
the rate of baryon outflow per unit solid angle Mq [g/s], and it is useful to rewrite T^ ec in terms of 
Mq. We substitute rii = rib = (Mq/ R 2 Tm p c) to equation (87), and then the decoupling condition 
recoil = R/c gives the maximum Lorentz factor of neutrons 



dec 



V 3 / \ V 3 , \ 1 / Q 



'—) - 300 (-i^- I ( -,. ) (91) 
r m p c J yi0 26 g s- 1 J \3x 10 6 cm/ 



4.2. Internal shocks 

The a-particles can also be destroyed later on, when internal shocks develop in the fireball. 
Lorentz factor r max fluctuates if the central engine is "noisy" during its operation < t < if, ~ 1 s. 
The fluctuations probably occur on timescales tq/c 10" 4 s and longer, up to if,. This leads to 
internal shocks (Rees Sz Meszaros 1994). Internal dissipation of the velocity fluctuations may give 
rise to the observed GRB, and this picture is plausible because it easily accounts for the observed 
short-timescale variability in GRB light curves. The amplitude of the fluctuations is described by 
the dimensionless rms of the Lorentz factor, A = Sr max /T max . At modest A ^ 1, the internal 
dissipation proceeds in a simple hierarchical manner (Beloborodov 2000). The internal collisions 
begin at 

*» - = 54 * loUA " ( w) 2 (stw) cm ' (95> 

where Ao [cm] is a minimum length scale of the fluctuations, which is likely comparable to r . The 
shortest fluctuations are dissipated first, and at R > Rq fluctuations with A > Ao are dissipated in 
the hierarchical order. 



4-2.1. n-ot collisions 

Suppose r max is small enough, so that the neutrons decouple from the ions after the acceleration 
stage and have same T n — r max . As soon as internal shocks develop in the ion component, the 
neutrons begin to drift with respect to the ions. Their relative velocity (3 pa A is sufficient for 
spallation reactions if A > 0.1 (and for A ~ 1 the pion production and multi-GeV neutrino emission 
takes place, see Meszaros & Rees 2000). 
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Let Lq = R nftm p c r max be the fireball kinetic luminosity per unit solid angle; then 

nb = 3l§F2 — • ( 96 ) 

rripC n i max 

The a-particles are destructed by neutrons with rate Y a « Y a Y n ribO- sp c and their lifetime in the 
fluid frame is t\h c = (Y^nb<7 sp c) _1 , 

rripC R r max 

Tlife = „ v t • ( 97 > 

It should be compared with the timescale of side expansion in the fluid frame, (-^/ cr max 

). if 

riife < {R/cT max ) at R = Rq then most of the a-particles are destroyed by the internal shocks. This 
condition reads 

V5 / r \ 1/5 /xx -1/5 



2Ao?77,pC 3 



4-2.2. a-ct collisions 

The a-particles acquire random energy (A 2 /2)4m p c 2 in internal shocks, which easily exceeds 
the spallation threshold. The condition riif e < (R/cT max ) for efficient a-a spallation is similar to 
equation (98) (with Y n replaced by Y a ). There is, however, one more condition. The shocked a- 
particles are cooled by Coulomb interactions with e~ (or e + ) on a timescale Tc ui> an d an efficient 
spallation requires T\tf e < tqouI i n addition to riif e < (R/cT ms , x ). 

The can be considered as targets at rest for the hot ions because their radiative losses 
(synchrotron and/or inverse Compton) are rapid compared to the expansion rate. In the fluid 
frame, hot ions with mass nii and charge Ze loose their random velocity f3 on timescale (Ginzburg 
k Syrovatskii 1964) 

Coul ^ ^ dr ^ 3Z 2 In Ao"T?T7. e cn e ' ^ ^ 

where n e = n_ + n + is the total density of e~ and e + , and In A ~ 20 is Coulomb logarithm. The 
internal shocks give j3 ~ A and, for a-particles (rrii = 4m p , Z = 2), we find riif e = (n a (j sp c) _1 < 
TCoul if 

< 2^3 rr^Osp ^ ^ 
n a 3mA m e cjt 

A crucial factor in this condition is the density. The postshock matter emits radiation, and 
pairs can be produced by 7-rays (7 + 7 —> e~ + e + ) that have energy hv > m e c 2 in the fluid frame. 
Suppose a fraction rj of the internal energy density {A 2 /2)nf,m p c 2 is converted into radiation and a 
fraction / of this radiation is above the threshold for pair production. A maximum pair density is 
evaluated assuming that all photons emitted above the threshold m e c 2 are converted into pairs, 

n b 2 m e 
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From equations (100) and (101) we conclude that the Coulomb losses of a-particles on e can 
prevent the efficient a-a spallation if fr\ > n = 3 x 10~ 3 AX a . In this case, only a fraction 
TCoui/Tiife ~ (i^/frj) °f the a-particles are destroyed. The released neutrons can then continue the 
spallation process via n-a collisions (§ 4.2.1). 

Equation (101) assumes an optical depth for 7-7 interactions r 77 > 1. We now check this 
assumption. The optical depth can be estimated « 0.1(u>i /m e c 2 )a T R/T 

max where w\ — 

fr)(A 2 /2)nbm p c 2 . This yields 

0.1/^^n 6 <7i.-^-. (102) 

Z illg 1 max 

It is easy to see that r 77 > 1 at all radii where efficient a-a spallation can take place [i.e. where 
riife < (R/cT max )} if fr) > 5 x lO" 4 . 



5. Conclusions 

1. Reactions of capture on nucleons, e~ + p — > n + v and e + + n — > p + &, operate in 
GRB central engines and set an equilibrium proton fraction Y e = n v j (n n + n p ). If the engine is an 
accretion disk around a black hole of mass M, the equilibrium Y e is established at accretion rates 
M > M cq w 1O 31 (q/O.1) 9 / 5 (M/M ) 6 / 5 g/s (eq. 48) where a = 0.01 - 0.1 is a viscosity parameter 
of the disk. 

2. Of great importance for the explosion dynamics is whether Y e < 0.5, i.e., there is an 
excess of neutrons. A general analysis shows that a neutrino-transparent matter has equilibrium 
Y e < 0.5 if ji > Q/2, and a similar condition for a z^-opaque matter reads [i > Q, where \i is the 
electron chemical potential in units of m e c 2 and Q = (m n — m p )/m e = 2.53. This condition is 
satisfied below a critical "neutronization" temperature T n (eqs. [14] and [23]). We find T < T n 
for plausible central engines of GRBs. In particular, accretion disks have a neutron excess at 
M > M n w 1O 31 (M/M ) 2 (a/O.l) g/s (eq. 50). 

3. Fireballs produced by neutron-rich engines should be also neutron rich. A major threat 
for neutrons in the escaping fireball is the neutrino flux from the central engine as absorption of a 
neutrino converts the neutron into a proton. This process is slower than the fireball expansion if 
the neutrino luminosity is below 10 53 erg/s. A neutrino luminosity above 10 53 erg/s would require 
a very powerful and i/-opaque central engine, which is possible. In this case, however, absorption 
of v by the fireball protons takes place as well as absorption of v by the neutrons. The balance 
between v and v absorptions establishes a new equilibrium Y e in the fireball, which depends on the 
emitted spectra of v and v and is likely below 0.5 (§ 2.3). 

4. As the fireball expands and cools, the ejected free nucleons tend to recombine into a- 
particles. This process competes, however, with rapid expansion and can freeze out. For this 
reason, nucleosynthesis is suppressed in fireballs with a high photon-to-baryon ratio 4> = n 7 /n{, 
(or, equivalently, high entropy per baryon s/k = 3.6^). We find that, in radial fireballs, more 
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than half of nucleons can recombine only if <j> < 3 x 10 4 (ro/3 x 10 6 ) - 85 where vq [cm] is the size 
of the central engine. In fireballs with parabolic collimation, the efficient recombination requires 
<fi < 6 x 10 5 (r /3 x 10 6 ) 09 . The typical GRB parameters <j) ~ 10 5 and r ~ 3 x 10 6 cm are just 
marginal for nucleosynthesis. 

5. Even in the case of efficient nucleon recombination, there are still leftover neutrons because 
of the neutron excess (Y e < 0.5). The minimum mass fraction of leftover neutrons is X n = 1 — 2Y e . 

6. The nucleosynthesis also produces deuterium, which is next abundant element after the free 
nucleons and ct-particles. Its typical mass fraction is a few per cent. The abundances of tritium, 
3 He, and all elements with mass number greater than 4 are negligible. 

7. Synthesized a-particles can be spalled later on, and then the populations of free neutrons 
and protons are increased. There are at least two possible mechanisms of spallation. (1) Energetic 
n-a collisions with a relative velocity (3 > f3 sp « 0.1 take place before the end of the fireball 
acceleration if <ft > 4 x 10 4 (L t h/I/) 3 / 4 . This mechanism works for fireballs with Lorentz factors 
Tmax > 160(/cTo/l MeV) (eq. 92) where To is an initial temperature (eq. 66). (2) Energetic n-a 
and a-a collisions occur when the fireball is reheated by internal shocks. This mechanism can 
be efficient at modest r m ax 

< 300(L n /10 52 ) 1/5 (eq. 98) where L n [erg/s] is the fireball kinetic 
luminosity per unit solid angle. 

The presence of a neutron component in GRB fireballs has quite spectacular implications. 
Beside making the fireball an interesting source of multi-GeV neutrinos (Derishev et al. 1999, 
Bahcall & Meszaros 2000, Meszaros k, Rees 2000) the neutrons survive and play a dramatic role 
for the explosion development at large radii R ~ 10 16 — 10 17 cm. When the fireball begins to 
decelerate as a result of the interaction with an external medium, neutrons continue to coast with 
a high Lorentz factor T n and form a leading front. One can easily show that this front has kinetic 
energy much larger (by a factor of X n T n T max ) than the rest-mass of the ambient medium. It leaves 
behind a relativistic trail loaded with the products of the neutron decay until the neutron front 
decays completely, which happens at R « 10 17 cm (about log X n r n r max « 10 times the mean decay 
radius). An external shock wave — the customary source of GRB afterglows — has to form in the 
neutron trail rather than in a normal ambient medium. The mechanism of the fireball deceleration 
in this situation is elaborated in Paper 2. 

The neutron features of GRB explosions appear inevitably in the standard fireball scenario. 
They would be absent, however, if the GRBs are produced by magnetized winds with extremely 
low baryon loading, where Poynting flux carries much more energy than matter (e.g. Usov 1994, 
Lyutikov & Blandford 2002). A detection or non-detection of neutron effects will constrain the 
level of baryon loading in GRBs. 

At the final stages of the preparation of this manuscript, a paper by Pruet, Woosley, & Hoffman 
(2002) appeared. They evaluate Y e for numerical accretion models of Popham et al. (1999). The 
results are consistent with our analysis in § 2 (note a slightly different definition of the viscosity 
parameter a OUT = [3/2]a Popham ). 
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